Code
#### STEP 0 - Loading the packages ####
library(sf) # simple features packages for handling vector GIS data
library(httr) # generic webservice package
library(tidyverse) # a suite of packages for data wrangling, transformation, plotting, ...
library(ows4R) # interface for OGC webservices
library(osmdata) # interface to access to Open street map data
library(dplyr)
library(rayrender) # 3D visualizations using R
library(elevatr) # interface to access elevation data
library(rayshader) # 3D visualizations using R
library(raster) # Package to manipulate rasters
library(lidR) # Package to manipulate lidar data
library(rgl) # Package to create 3D interactive graphics
library(rgdal) # Geospatial data
# Define central point
central_point = c(5.664740537573191,51.98513909689059)
central_point|>
st_point() |>
st_sfc(crs = 4326) |>
st_transform(st_crs("EPSG:28992")) ->
central_point
# Define area of interest (Buffer around central point)
aoi_distance = units::as_units(400,"m")
aoi_buffer = st_buffer(central_point,aoi_distance)
aoi_bbox <- st_bbox(aoi_buffer)
bbox_string = paste(aoi_bbox[1],aoi_bbox[2],aoi_bbox[3],aoi_bbox[4],sep = ",")
# Download 3dBag from the area:
wfs_3dbag <- "https://data.3dbag.nl/api/BAG3D_v2/wfs"
wfs_3dbag_client <- WFSClient$new(wfs_3dbag,
serviceVersion = "2.0.0")
url <- parse_url(wfs_3dbag)
url$query <- list(service = "wfs",
version = "2.0.0",
request = "GetFeature",
typename = "bag_tiles_3k",
bbox = bbox_string,
srsName = "EPSG:28992",
outputFormat = "gml32"
)
# Get tiles to download
request <- build_url(url)
tiles_3d_aoi <- read_sf(request)
# Download files
folder_download <- "data_nl_3dbag"
dir.create(folder_download)
# Downloads all tiles in the area
buildings <- lapply(tiles_3d_aoi$tile_id, function(x){
file_name <- paste0("3dbag_v210908_fd2cee53_",x,".gpkg")
url_download <- paste0("https://data.3dbag.nl/gpkg/v210908_fd2cee53/",file_name)
GET(url_download, write_disk(file.path(folder_download,file_name),overwrite = TRUE))
st_read(file.path(folder_download,file_name),layer = "lod22_3d", quiet=TRUE)
}
)
# Merges all buildings into a single element
# This only works because all files are using the same Coordinate reference system !!!
city_buildings <- do.call(rbind,buildings)
# Change the name of the geometry [Dutch to English... (sorry)]
st_geometry(city_buildings)<-"geometry"
city_buildings|> st_transform(st_crs(central_point)) -> city_buildings
# Select only the buildings within the area of interests
city_buildings |>
# Defines centroid of each bulding
mutate(centroid = st_centroid(geometry)) |>
# Computes distance to central point
mutate(dist = st_distance(centroid, central_point)) |>
# Selects buildings within the AIO
filter(dist < aoi_distance) -> city_buildings_aoi
# Loading the elevation raster
elevation_data = elevatr::get_elev_raster(
locations = st_buffer(st_centroid(st_union(city_buildings_aoi)),
dist=aoi_distance),z=14)
#Get the bounding box for the scene
scene_bbox = st_bbox(city_buildings_aoi)
#Crop the elevation data to that bounding box
cropped_data = raster::crop(elevation_data, scene_bbox)
# Create a matrix from the elevation
elevation_matrix = raster_to_matrix(cropped_data)
# Render
elevation_matrix %>%
height_shade() %>%
plot_3d(elevation_matrix,baseshape = "rectangle")
render_multipolygonz(city_buildings_aoi,
extent = raster::extent(cropped_data),
baseshape = "rectangle",
clear_previous = TRUE,color = "grey",offset = 0,zscale = 1,
heightmap = elevation_matrix)
rglwidget()3D plot

📍 (Gaia 3rd floor)
📍 (Gaia 3rd floor)